Project background

Essentially all individuals experience social stress, which may be acute or chronic. Individuals differ in their behavioral and physiological responses to stressors, but what factors cause these individual differences is less clear. Manipulations to the maternal microbiome and maternal stress can have sex-specific effects on offspring social behaviors and may help explain variation in behavioral responses to prolonged stress. My research project investigates how maternal environment affects differences in offspring social behavior during long bouts of stress in Siberian hamsters (Phodopus sungorus, figure 1).

**Figure 1.** Image of Siberian hamster *Phodopus sungorus*. Image attribution: Salzgeber, P. (2007). Phodopus sungorus. In Wikimedia Commons. https://commons.wikimedia.org/wiki/File:PhodopusSungorus_1.jpg

Figure 1. Image of Siberian hamster Phodopus sungorus. Image attribution: Salzgeber, P. (2007). Phodopus sungorus. In Wikimedia Commons. https://commons.wikimedia.org/wiki/File:PhodopusSungorus_1.jpg

Research questions

Question: How does early development (i.e., maternal environment) affect whether offspring vary in their social behaviors during the course of a long-bout social stressor?
Hypothesis: Maternal stress during offspring development has sex-specific effects on consistency in individual behavior during a long-bout social stressor
Predictions:
1. Female offspring from Stress Only mothers will display increased aggression throughout a long-bout stressor.
2. Other offspring will not display consistent behaviors throughout the long-bout stressor.

Methods

Maternal treatment

Pregnant females (n=34) were exposed to one of four treatments for ten days during pregnancy: control, antibiotic only, stress only, antibiotic + stress (figure 2)
**Figure 2.** Maternal treatment graphic depicting the four treatment groups.

Figure 2. Maternal treatment graphic depicting the four treatment groups.


Offspring behavioral trials

Juvenile offspring (n=14 males, 11 females) were exposed to a 15 minute social stressor following a resident-intruder paradigm. The intruder was placed into the focal individual’s (i.e., offspring’s) cage for 15 minutes, and offspring behaviors towards the intruder were scored (figure 3).
The first 5 minutes (t0-t5) and the last five minutes (t10-t15) of behaviors were scored and analyzed.
**Figure 3.**Offspring in late adolescence (51-56 days old) were exposed to a resident-intruder trial as above. The intruder was placed in the focal individual’s cage for 15 minutes. Time scored is shown by red line.

Figure 3.Offspring in late adolescence (51-56 days old) were exposed to a resident-intruder trial as above. The intruder was placed in the focal individual’s cage for 15 minutes. Time scored is shown by red line.


Four behaviors were scored: attack, investigation, escape, and grooming (figure 4).

**Figure 4.** Actual screenshots of video data depicting the four social behaviors of interest.

Figure 4. Actual screenshots of video data depicting the four social behaviors of interest.

Data cleaning

Below is the code used to clean the and prepare the data for analyses. Several columns were removed from the data frame because they were not used in this project (these were columns imported into the excel sheet from the original project in 2019).

# Load data, clean names, view structure ####
hammies <- read_csv("../final_project/IU106_Social Behavior First and Last Five Minutes_Formatted_6-6-23.csv")
hammies <- janitor::clean_names(hammies)
names(hammies)
str(hammies)

# Remove columns that were not used in this project ####
columns_removed <- c("id", "aggression_number_assigned", "comp1", "comp2", "comp3", 
                     "aggression_score_freq", "aggression_score_duration")
hammies <- dplyr::select(hammies, -columns_removed)

# Create trial timing column, set appropriate columns as factors ####
tidy_hammies <- hammies %>% 
  mutate(trial_timing = as.factor(first_or_last_five_minutes)) %>% 
  mutate(treatment = as.factor(maternal_treatment)) %>% 
  mutate(juvenile_id = as.factor(juvenile_id)) %>% 
  dplyr::select(-first_or_last_five_minutes)

# Fix spelling of all "occurence" columns ####
tidy_hammies <- tidy_hammies %>% 
  mutate(attack_number_of_occurrences = attack_number_of_occurences) %>%
  mutate(chase_number_of_occurrences = chase_number_of_occurences) %>%
  mutate(groom_number_of_occurrences = groom_number_of_occurences) %>%
  mutate(investigation_number_of_occurrences = investigation_number_of_occurences) %>%
  mutate(jump_number_of_occurrences = jump_number_of_occurences) %>% 
  dplyr::select(-contains("occurences"))

str(tidy_hammies)

Modeling Juvenile Behavior

We analyzed the number (i.e., frequency) and total duration of behavioral bouts in the first and last five minutes (t0-t5 vs. t10-t15).
For each model, we first tested for a three-way interaction between maternal treatment, offspring sex, and trial timing. If no significant three-way interaction was detected, we tested for a two-way interaction between offspring sex and maternal treatment. Only significant interactions were retained in the model; non-significant interaction terms were removed.

Summary of models

All models had two random effects: intruder identity and resident (offspring) identity, to control for individual differences in behavior. All frequency models (GLMM’s) used the bobyqa optimizer from the lme4 package for model fitting.

Table 1. Summary of models run for each response variable.
Response Variable Distribution Fixed Effects
Attack Frequency Poisson, Log link function Offspring sex, maternal treatment, trial timing, SI- CORT, offspring weight, offspring sex*maternal treatment interaction
Attack Duration Gaussian, Identity link function Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Investigation Frequency Poisson, Log link function Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sexmaternal treatmenttrial timing
Investigation Duration Gaussian, Identity link function Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Escape Frequency Poisson, Log link function Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sexmaternal treatmenttrial timing
Escape Duration Gaussian, Identity link function Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Grooming Frequency Poisson, Log link function Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Grooming Duration Gaussian, Identity link function Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sex* maternal treatment

Agression frequency

Results

Our first model analyzes the effects of maternal treatment, offspring sex, and video timing on the frequency of aggressive behaviors. There was no significant three-way interaction, but we did detect a significant interaction between maternal treatment and offspring sex. Below is the code for model A.1.

modelA.1 <- tidy_hammies %>% 
  glmer(attack_number_of_occurrences ~ sex*treatment+trial_timing+scale(si_cort)
        +weight_at_trial+(1|intruder_id)+(1|juvenile_id),family=poisson,
        glmerControl(optimizer="bobyqa"),data=.)
broom.mixed::tidy(modelA.1) %>% 
  kable("pipe", digits = 3, caption = "**Table 2.** Summary of model A.1, analyzing the interactive effects of maternal treatment and offspring sex on frequency of attack behaviors in the first and last five minutes of the resident-intruder trial.") %>% 
  kable_styling() %>%
  kableExtra::row_spec(.,c(4:6, 10), bold=TRUE) %>% 
  kableExtra::column_spec(column = 1:7, width = "200px")
Table 2. Summary of model A.1, analyzing the interactive effects of maternal treatment and offspring sex on frequency of attack behaviors in the first and last five minutes of the resident-intruder trial.
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 0.589 0.706 0.834 0.404
fixed NA sexMale 0.430 0.373 1.153 0.249
fixed NA treatmentAntibiotic Only 0.604 0.388 1.555 0.120
fixed NA treatmentControl -2.299 1.050 -2.190 0.029
fixed NA treatmentStress Only 0.852 0.384 2.221 0.026
fixed NA trial_timingLast -0.321 0.153 -2.103 0.035
fixed NA scale(si_cort) 0.023 0.101 0.228 0.819
fixed NA weight_at_trial 0.014 0.022 0.645 0.519
fixed NA sexMale:treatmentAntibiotic Only -0.551 0.485 -1.137 0.256
fixed NA sexMale:treatmentControl 2.323 1.106 2.101 0.036
fixed NA sexMale:treatmentStress Only -0.784 0.482 -1.626 0.104
ran_pars juvenile_id sd__(Intercept) 0.145 NA NA NA
ran_pars intruder_id sd__(Intercept) 0.000 NA NA NA

Model A.1 is a poisson mixed model fitted to predict number of attack bouts. 

tidy_hammies %>% 
  ggplot(aes(x=interaction(trial_timing, sex),
             y = attack_number_of_occurrences))+
  geom_boxplot(aes(fill=treatment))+
  labs(x= "Trial timing", 
       y = "Number of attacks", 
       fill = "Maternal treatment")+
  theme(axis.line = element_line(color="black"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank())
**Figure 5.** Plot of aggression frequency by sex and trial timing.

Figure 5. Plot of aggression frequency by sex and trial timing.


We then used the lsmeans function from the emmeans R package to conduct pairwise comparisons of the maternal treatment and sex interactions. lsmeans uses a corrected p-value since to minimize the likelihood of a type I error, since we are performing a large number of comparisons. After running the lsmeans, none of our pairwise comparisons are statistically significant using a p-value of 0.05.

#pairwise comparisons with corrected p-value
lsmeans(modelA.1, pairwise~treatment, type = "response") # because there was a significant 2-way interaction, we are more interested in the interaction below.
lsmeans(modelA.1, pairwise~treatment*sex, type = "response") 
lsmeans_output <- lsmeans(modelA.1, pairwise~treatment*sex, type = "response")
lsmeans_contrasts <- lsmeans_output$contrasts
lsmeans_contrasts_df <- as.data.frame(lsmeans_contrasts)
lsmeans_contrasts_df %>% 
  kable("pipe", digits = 3, caption = "**Table 3.** lsmeans comparisons for attack frequency.") %>% 
  kable_styling() %>%
  kableExtra::row_spec(.,14, bold=TRUE) %>% 
  kableExtra::column_spec(column = 1:7, width = "200px")
Table 3. lsmeans comparisons for attack frequency.
contrast estimate SE df null z.ratio p.value
(Antibiotic + Stress Female) / Antibiotic Only Female 0.547 0.212 Inf 1 -1.555 0.777
(Antibiotic + Stress Female) / Control Female 9.964 10.463 Inf 1 2.190 0.358
(Antibiotic + Stress Female) / Stress Only Female 0.427 0.164 Inf 1 -2.221 0.339
(Antibiotic + Stress Female) / (Antibiotic + Stress Male) 0.651 0.243 Inf 1 -1.153 0.945
(Antibiotic + Stress Female) / Antibiotic Only Male 0.617 0.217 Inf 1 -1.370 0.871
(Antibiotic + Stress Female) / Control Male 0.635 0.219 Inf 1 -1.316 0.893
(Antibiotic + Stress Female) / Stress Only Male 0.608 0.231 Inf 1 -1.307 0.897
Antibiotic Only Female / Control Female 18.225 19.168 Inf 1 2.760 0.105
Antibiotic Only Female / Stress Only Female 0.780 0.273 Inf 1 -0.710 0.997
Antibiotic Only Female / (Antibiotic + Stress Male) 1.190 0.445 Inf 1 0.465 1.000
Antibiotic Only Female / Antibiotic Only Male 1.129 0.389 Inf 1 0.352 1.000
Antibiotic Only Female / Control Male 1.161 0.395 Inf 1 0.440 1.000
Antibiotic Only Female / Stress Only Male 1.112 0.418 Inf 1 0.282 1.000
Control Female / Stress Only Female 0.043 0.045 Inf 1 -2.975 0.059
Control Female / (Antibiotic + Stress Male) 0.065 0.067 Inf 1 -2.648 0.139
Control Female / Antibiotic Only Male 0.062 0.064 Inf 1 -2.709 0.120
Control Female / Control Male 0.064 0.066 Inf 1 -2.672 0.131
Control Female / Stress Only Male 0.061 0.063 Inf 1 -2.709 0.120
Stress Only Female / (Antibiotic + Stress Male) 1.525 0.588 Inf 1 1.095 0.958
Stress Only Female / Antibiotic Only Male 1.447 0.485 Inf 1 1.103 0.956
Stress Only Female / Control Male 1.489 0.513 Inf 1 1.155 0.944
Stress Only Female / Stress Only Male 1.425 0.529 Inf 1 0.955 0.980
(Antibiotic + Stress Male) / Antibiotic Only Male 0.949 0.285 Inf 1 -0.176 1.000
(Antibiotic + Stress Male) / Control Male 0.976 0.308 Inf 1 -0.077 1.000
(Antibiotic + Stress Male) / Stress Only Male 0.935 0.299 Inf 1 -0.212 1.000
Antibiotic Only Male / Control Male 1.029 0.307 Inf 1 0.096 1.000
Antibiotic Only Male / Stress Only Male 0.985 0.284 Inf 1 -0.052 1.000
Control Male / Stress Only Male 0.958 0.317 Inf 1 -0.131 1.000

Model assumptions and validation

Below is the code my PI initially used to check model assumptions and check for overdispersion.

# * * Checking model assumptions ------------------------------------------
qqnorm(residuals(modelA.1))
qqline(residuals(modelA.1))
shapiro.test(residuals(modelA.1))


# * * Check overdispersion ------------------------------------------------
overdisp_fun <- function(modelA.1) {
  rdf <- df.residual(modelA.1)
  rp <- residuals(modelA.1,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun
overdisp_fun(modelA.1)

sum(residuals(modelA.1,type="pearson")^2)/(nrow(tidy_hammies)-length(fixef(modelA.1)-1))

vif(modelA.1)

Using the performance package we learned about in class, I simplified the overdispersion code, and used check_model to visualize vif and overdispersion.

check_model(modelA.1)
**Figure 6.** Output of performance::check_model function for model A.1.

Figure 6. Output of performance::check_model function for model A.1. 

performance::check_overdispersion(modelA.1)
## # Overdispersion test
## 
##        dispersion ratio =  1.341
##   Pearson's Chi-Squared = 46.932
##                 p-value =  0.086

Aggression Duration

We fitted a linear mixed model for the duration of attack behaviors in the first and last five minutes of the resident-intruder trial. The data were square-root transformed to assume a normal distribution. There were no significant interactions to predict attack behavior duration.

modelA.2 <- lmer(sqrt(attack_total_duration)~ sex+treatment+trial_timing+
                   scale(si_cort)+weight_at_trial+(1|intruder_id)+(1|juvenile_id),
                 data = tidy_hammies)
#summary(modelA.2)
broom.mixed::tidy(modelA.2) %>% 
  kable("pipe", digits = 3) %>% 
  kable_styling() %>%
  kableExtra::column_spec(column = 1:7, width = "200px")
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 3.608 1.390 2.597 17.398 0.019
fixed NA sexMale 0.749 0.452 1.656 17.000 0.116
fixed NA treatmentAntibiotic Only 0.324 0.600 0.540 17.000 0.596
fixed NA treatmentControl -0.256 0.586 -0.437 17.000 0.668
fixed NA treatmentStress Only 0.397 0.631 0.630 17.000 0.537
fixed NA trial_timingLast -0.009 0.299 -0.031 23.000 0.976
fixed NA scale(si_cort) 0.180 0.251 0.717 17.000 0.483
fixed NA weight_at_trial -0.069 0.043 -1.584 17.000 0.132
ran_pars juvenile_id sd__(Intercept) 0.677 NA NA NA NA
ran_pars intruder_id sd__(Intercept) 0.000 NA NA NA NA
ran_pars Residual sd__Observation 1.034 NA NA NA NA

Summary and Conclusions

  • Overall, both frequency and duration of social behaviors decrease over a 15 minute trial
  • Maternal treatment has significant effects on attack frequency.
  • Maternal treatment has sex-specific effects on investigation frequency and duration, escape frequency, and grooming duration.
  • Trial timing has significant effects on escape duration.
  • This research suggests that early maternal experience affects behavior during long bouts of stress, and effects can be sex-specific
  • Mechanisms driving changes in social behaviors during long bouts of stress need to be explored further
modelC.3 <- lmer(jump_total_duration_s ~sex+treatment+trial_timing+scale(si_cort)+
                   weight_at_trial+(1|intruder_id)+(1|juvenile_id), 
                 data = tidy_hammies)
# relevel treatment to control
tidy_hammies <- tidy_hammies %>% 
  mutate(treatment = relevel(treatment, ref = "Control"))
tidy_hammies %>% 
  ggplot(aes(x=trial_timing, 
             y = jump_total_duration_s))+
  geom_boxplot(aes(fill = sex))+
  labs(x= "Trial timing", 
       y = "Duration of jump behaviors", 
       fill = "Sex")+
  theme(axis.line = element_line(color="black"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank())
**Figure 7.** Plot of escape duration by sex and trial timing. This graph shows a huge increase in the duration of escape behaviors in the last five minutes of the trial across maternal treatment groups and across sex.

Figure 7. Plot of escape duration by sex and trial timing. This graph shows a huge increase in the duration of escape behaviors in the last five minutes of the trial across maternal treatment groups and across sex.